# ------------------------------------------------------------
# MIND THE GAP:  Replication Script.  
# Date: 2025-10-29
# ------------------------------------------------------------
# REPLICATION INSTRUCTIONS:
# -------------------------
# 1. SET YOUR WORKING DIRECTORY
# 2. PLACE YOUR DATA FILE in that directory
# 3. UPDATE THE DATA FILENAME if different
# 4. RUN THE ENTIRE SCRIPT (Ctrl+Shift+Enter in RStudio)
#
#
# Expected runtime: 1-2 minutes
# Required R version: >= 4.0.0
# ============================================================

# ============================================================
# USER CONFIGURATION - EDIT THESE LINES BEFORE RUNNING
# ============================================================

# Set your working directory (where data file is located)
setwd(getwd())  # Default: current directory. Change to your path: setwd("C:/Users/YourName/Documents")

# Set your data filename
#data_filename <- "NO26467145_Power_of_Party_data_Final (1).sav"  # REMOVE # AND CHANGE IF NEEDED

# Verify setup
if (!file.exists(data_filename)) {
  stop("\n*** ERROR: Data file not found! ***\n",
       "Looking for: ", file.path(getwd(), data_filename), "\n",
       "Please check:\n",
       "  1. Working directory is correct \n",
       "  2. Filename matches exactly \n",
       "  3. File is in the specified directory \n")
}
# ============================================================

# Sections:
# I. Data import & cleaning
# II. Descriptives (console)
# III. Hypothesis tests H1-H3 + omnibus + H5-H6 (console)
# IV. Exploratory - tests (console)
# V. Post hoc power & MDES (console)
# VI. Appendix: Dark personality (H4, H7) tables
# VII. Robustness with covariates (R0-R6)
# Z. Export Word file with all tables
# F. Main figures (H1-H3)

# --- clean start ---
rm(list = ls()); gc()

# --- end clean start ---

# encoding
try({ if (.Platform$OS.type == "windows") Sys.setlocale("LC_CTYPE", "English_United States.utf8") }, silent = TRUE)

# ---- Packages: install if missing, then load (quietly) ----
need <- c(
  "dplyr","haven","car","effectsize","pwr","broom","ggplot2",
  "mgcv","emmeans","patchwork","tibble","officer","flextable"
)
to_install <- need[!sapply(need, requireNamespace, quietly = TRUE)]
if (length(to_install)) install.packages(to_install)
invisible(lapply(need, function(pk) suppressPackageStartupMessages(library(pk, character.only = TRUE))))


# Using R defaults (treatment contrasts) for coefficient interpretation
options(contrasts = c("contr.treatment", "contr.poly"))

# ============================================================
# STEP I. Data import, variable construction, cleaning
# ============================================================
dat <- read_sav(data_filename) #NO26467145_Power_of_Party_data_Final (1).sav


hypotheses0 <- dat %>%
  mutate(
    response = coalesce(B1, B2, B3, B4, B5, B6),
    disruption_severity = case_when(
      !is.na(B1) | !is.na(B3) | !is.na(B5) ~ "Mild",
      !is.na(B2) | !is.na(B4) | !is.na(B6) ~ "Severe",
      TRUE ~ NA_character_
    ),
    party_condition = case_when(
      !is.na(B1) | !is.na(B2) ~ "No Party",
      !is.na(B3) | !is.na(B4) ~ "Outparty",
      !is.na(B5) | !is.na(B6) ~ "Inparty",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(response))

hypotheses0 <- hypotheses0 %>%
  mutate(
    partisanship = case_when(
      (pid3 %in% c(1, 2) & A2 == 1) ~ "Strong",
      (pid3 %in% c(1, 2) & A2 == 2) ~ "Moderate",
      (pid3 == 3 & A3 %in% c(1, 2)) ~ "Moderate",
      TRUE ~ NA_character_
    ),
    party_family = case_when(
      pid3 == 1 ~ "Democrat",
      pid3 == 2 ~ "Republican",
      pid3 == 3 & A3 == 1 ~ "Democrat",
      pid3 == 3 & A3 == 2 ~ "Republican",
      TRUE ~ NA_character_
    ),
    partisanship_dummy = ifelse(partisanship == "Strong", 1, ifelse(partisanship == "Moderate", 0, NA))
  )

hypotheses_partisans <- hypotheses0 %>%
  filter(!is.na(partisanship), QAttention == 4) %>%
  mutate(
    party_condition = factor(party_condition, levels = c("No Party","Inparty","Outparty")),
    disruption_severity = factor(disruption_severity, levels = c("Mild","Severe")),
    party_family = factor(party_family, levels = c("Democrat","Republican"))
  )

# ============================================================
# STEP II. Descriptives - console prints
# ============================================================

cat("\n--- Descriptives ---\n")

N_all <- nrow(hypotheses0)
N_inattentive <- sum(hypotheses0$QAttention != 4 & !is.na(hypotheses0$QAttention))
N_partisans <- sum(!is.na(hypotheses0$partisanship))
N_nonpartisans <- N_all - N_partisans
N_retained <- nrow(hypotheses_partisans)

cat(sprintf("Total with non-NA response: %d\n", N_all))
cat(sprintf("Inattentive (QAttention != 4): %d\n", N_inattentive))
cat(sprintf("Non-partisans (incl. true independents): %d\n", N_nonpartisans))
cat(sprintf("Retained (attentive partisans): %d\n", N_retained))

tab_npct <- function(x) {
  x_no_na <- x[!is.na(x)]
  df <- as.data.frame(table(x_no_na), stringsAsFactors = FALSE)
  names(df) <- c("Category","N")
  df <- df %>% mutate(Percent = round(100 * N / sum(N), 1))
  df[order(-df$N), ]
}

cat("\nPartisanship (Strong vs Moderate) among retained:\n")
print(tab_npct(hypotheses_partisans$partisanship))

cat("\nParty family among retained:\n")
print(tab_npct(hypotheses_partisans$party_family))

cat("\nGender distribution among retained:\n")
gender_vec <- if (haven::is.labelled(hypotheses_partisans$Q_Gender)) {
  haven::as_factor(hypotheses_partisans$Q_Gender, levels = "labels")
} else hypotheses_partisans$Q_Gender
print(tab_npct(gender_vec))

cat("\nAge summary among retained:\n")
print(summary(hypotheses_partisans$age))
cat("\nAge (10-year bins) among retained:\n")
age_bins <- cut(hypotheses_partisans$age, breaks = c(18,25,35,45,55,65,75,85,Inf), right = FALSE)
print(tab_npct(age_bins))

income_clean <- if (haven::is.labelled(hypotheses_partisans$profile_gross_household)) {
  haven::as_factor(hypotheses_partisans$profile_gross_household, levels = "labels")
} else hypotheses_partisans$profile_gross_household
if (is.factor(income_clean)) {
  income_clean <- droplevels(income_clean[income_clean != "Don't know/Prefer not to say"])
} else {
  income_clean[income_clean == 99] <- NA
}
cat("\nIncome distribution among retained:\n")
print(tab_npct(income_clean))

race_vec <- if (haven::is.labelled(hypotheses_partisans$race)) {
  haven::as_factor(hypotheses_partisans$race, levels = "labels")
} else hypotheses_partisans$race
cat("\nRace/ethnicity among retained:\n")
print(tab_npct(race_vec))

cat("\nCell counts by party_condition:\n")
print(hypotheses_partisans %>% count(party_condition, name = "N"))
cat("\nCell counts by disruption_severity:\n")
print(hypotheses_partisans %>% count(disruption_severity, name = "N"))
cat("\nCell counts by party_condition - disruption_severity:\n")
print(hypotheses_partisans %>% count(party_condition, disruption_severity, name = "N") %>% arrange(party_condition, disruption_severity))

# ============================================================
# STEP III. Hypotheses H1-H3, omnibus, H5-H6 - console prints
# ============================================================

cat("\n--- Hypothesis testing (attentive partisans only) ---\n")

hypotheses_partisans <- hypotheses_partisans %>%
  mutate(
    party_family = case_when(
      pid3 == 1 ~ "Democrat",
      pid3 == 2 ~ "Republican",
      pid3 == 3 & A3 == 1 ~ "Democrat",
      pid3 == 3 & A3 == 2 ~ "Republican",
      TRUE ~ NA_character_
    )
  )

no_party <- filter(hypotheses_partisans, party_condition == "No Party")
inparty  <- filter(hypotheses_partisans, party_condition == "Inparty")
outparty <- filter(hypotheses_partisans, party_condition == "Outparty")

t_H1_out_vs_no  <- t.test(outparty$response, no_party$response, alternative = "two.sided")
t_H1_in_vs_no   <- t.test(inparty$response,  no_party$response,  alternative = "two.sided")
t_H2            <- t.test(outparty$response, inparty$response, alternative = "greater")

mild   <- filter(hypotheses_partisans, disruption_severity == "Mild")
severe <- filter(hypotheses_partisans, disruption_severity == "Severe")
t_H3_all <- t.test(severe$response, mild$response, alternative = "greater")
t_H3_no  <- with(no_party, t.test(response[disruption_severity=="Severe"],
                                  response[disruption_severity=="Mild"], alternative="greater"))
t_H3_out <- with(outparty, t.test(response[disruption_severity=="Severe"],
                                  response[disruption_severity=="Mild"], alternative="greater"))
t_H3_in  <- with(inparty,  t.test(response[disruption_severity=="Severe"],
                                  response[disruption_severity=="Mild"], alternative="greater"))

cat("\nH1 t-test: Outparty vs No Party (two-sided)\n"); print(t_H1_out_vs_no)
cat("\nH1 t-test: Inparty vs No Party (two-sided)\n");  print(t_H1_in_vs_no)
cat("\nH2 t-test: Outparty > Inparty (one-sided)\n");   print(t_H2)
cat("\nH3 t-test: Severe > Mild (one-sided)\n");        print(t_H3_all)
cat("\nH3 within-group: No Party (Severe > Mild)\n");   print(t_H3_no)
cat("\nH3 within-group: Outparty (Severe > Mild)\n");   print(t_H3_out)
cat("\nH3 within-group: Inparty (Severe > Mild)\n");    print(t_H3_in)

a_H1 <- aov(response ~ party_condition, data = hypotheses_partisans)
cat("\nH1 omnibus ANOVA (equal variances):\n"); print(summary(a_H1))

w_H1 <- oneway.test(response ~ party_condition, data = hypotheses_partisans, var.equal = FALSE)
cat("\nH1 omnibus ANOVA (Welch):\n"); print(w_H1)

a_H2 <- aov(response ~ factor(party_condition, levels = c("Inparty","Outparty")),
            data = filter(hypotheses_partisans, party_condition %in% c("Inparty","Outparty")))
cat("\nH2 omnibus ANOVA (2 levels: Inparty vs Outparty):\n"); print(summary(a_H2))

a_H3 <- aov(response ~ disruption_severity, data = hypotheses_partisans)
cat("\nH3 omnibus ANOVA (2 levels: Mild vs Severe):\n"); print(summary(a_H3))

op_contr <- options(contrasts = c("contr.sum","contr.poly"))
m_2x3 <- lm(response ~ party_condition * disruption_severity, data = hypotheses_partisans)
a3_2x3 <- car::Anova(m_2x3, type = "III")
options(op_contr)
cat("\nFull 2x3 factorial (Type III ANOVA):\n"); print(a3_2x3)

# ============================================================
# STEP IV. Explorative - proportion tests - console prints
# ============================================================

cat("\n--- Explorative expansion: explicit condonance (x^2 one-tailed) ---\n")

hypotheses_partisans <- hypotheses_partisans %>%
  mutate(
    party_family = dplyr::case_when(
      pid3 == 1 ~ "Democrat",
      pid3 == 2 ~ "Republican",
      pid3 == 3 & A3 == 1 ~ "Democrat",
      pid3 == 3 & A3 == 2 ~ "Republican",
      TRUE ~ NA_character_
    ),
    condone = as.numeric(response) >= 5
  )

fmt_p <- function(p){
  if (is.na(p)) return(NA_character_)
  if (p >= 1e-4) sprintf("%.6f", p) else format(p, digits = 16, scientific = TRUE)
}
scrub_inf <- function(x){
  x[is.infinite(x) & x < 0] <- -1
  x[is.infinite(x) & x > 0] <-  1
  x
}
prop_diff_stats <- function(df, sel1, sel2, alt = "greater") {
  x <- c(sum(df$condone[sel1], na.rm = TRUE), sum(df$condone[sel2], na.rm = TRUE))
  n <- c(sum(sel1, na.rm = TRUE),            sum(sel2, na.rm = TRUE))
  if (any(n < 1)) return(list(est = NA_real_, p = NA_real_, ci = c(NA_real_, NA_real_), n = n))
  tst <- suppressWarnings(prop.test(x = x, n = n, alternative = alt, correct = FALSE))
  est <- unname(tst$estimate[1] - tst$estimate[2])
  ci  <- scrub_inf(unname(tst$conf.int))
  list(est = est, p = tst$p.value, ci = ci, n = n)
}

G <- list(
  "All Partisans"          = function(d) rep(TRUE, nrow(d)),
  "Strong Partisans"       = function(d) d$partisanship == "Strong",
  "Moderate Partisans"     = function(d) d$partisanship == "Moderate",
  "All Republicans"        = function(d) d$party_family == "Republican",
  "Strong Republicans"     = function(d) d$party_family == "Republican" & d$partisanship == "Strong",
  "Moderate Republicans"   = function(d) d$party_family == "Republican" & d$partisanship == "Moderate",
  "All Democrats"          = function(d) d$party_family == "Democrat",
  "Strong Democrats"       = function(d) d$party_family == "Democrat" & d$partisanship == "Strong",
  "Moderate Democrats"     = function(d) d$party_family == "Democrat" & d$partisanship == "Moderate"
)

cat("\n[H2] One-tailed chi-square: Outparty - Inparty\n")
for (lbl in names(G)) {
  dd <- hypotheses_partisans[G[[lbl]](hypotheses_partisans), , drop = FALSE]
  st <- prop_diff_stats(dd, dd$party_condition == "Outparty", dd$party_condition == "Inparty", alt = "greater")
  cat(sprintf("%-22s  est = %.4f, p = %s, CI = [%.4f, %.4f], n1=%d, n2=%d\n",
              lbl, st$est, fmt_p(st$p), st$ci[1], st$ci[2], st$n[1], st$n[2]))
}

cat("\n[H1] One-tailed chi-square with No Party first\n")
for (lbl in names(G)) {
  dd <- hypotheses_partisans[G[[lbl]](hypotheses_partisans), , drop = FALSE]
  st1 <- prop_diff_stats(dd, dd$party_condition == "No Party", dd$party_condition == "Outparty", alt = "less")
  st2 <- prop_diff_stats(dd, dd$party_condition == "No Party", dd$party_condition == "Inparty",  alt = "greater")
  cat(sprintf("%-22s  NoParty-Out: est=%.4f, p=%s, CI=[%.4f, %.4f];  NoParty-In: est=%.4f, p=%s, CI=[%.4f, %.4f]  (nNo=%d / nOut=%d / nIn=%d)\n",
              lbl,
              st1$est, fmt_p(st1$p), st1$ci[1], st1$ci[2],
              st2$est, fmt_p(st2$p), st2$ci[1], st2$ci[2],
              sum(dd$party_condition=="No Party"),
              sum(dd$party_condition=="Outparty"),
              sum(dd$party_condition=="Inparty")))
}

cat("\n[H3] One-tailed chi-square: Severe - Mild\n")
for (lbl in names(G)) {
  dd <- hypotheses_partisans[G[[lbl]](hypotheses_partisans), , drop = FALSE]
  st <- prop_diff_stats(dd, dd$disruption_severity == "Severe", dd$disruption_severity == "Mild", alt = "greater")
  cat(sprintf("%-22s  est = %.4f, p = %s, CI = [%.4f, %.4f], nSev=%d, nMild=%d\n",
              lbl, st$est, fmt_p(st$p), st$ci[1], st$ci[2], st$n[1], st$n[2]))
}

# ============================================================
# STEP V. Post hoc power and MDES - console prints
# ============================================================

cat("\n--- Post hoc power and MDES ---\n")

alpha <- 0.05
harmonic_mean <- function(n_vec) { n_vec <- n_vec[n_vec > 0]; length(n_vec) / sum(1 / n_vec) }

power_oneway <- function(data, y, group, label) {
  dsub <- data %>% filter(!is.na(.data[[y]]), !is.na(.data[[group]]))
  dsub[[group]] <- droplevels(as.factor(dsub[[group]]))
  k <- nlevels(dsub[[group]])
  ns <- table(dsub[[group]])
  n_h <- harmonic_mean(as.numeric(ns))
  fit <- aov(reformulate(group, response = y), data = dsub)
  tab <- summary(fit)[[1]]
  SS_eff <- as.numeric(tab[1, "Sum Sq"])
  SS_res <- as.numeric(tab[nrow(tab), "Sum Sq"])
  eta2   <- SS_eff / (SS_eff + SS_res)
  f_obs  <- sqrt(eta2 / (1 - eta2))
  pw       <- pwr.anova.test(k = k, n = n_h, f = f_obs, sig.level = alpha)$power
  mdes_obj <- pwr.anova.test(k = k, n = n_h, sig.level = alpha, power = 0.80)
  f_mdes   <- as.numeric(mdes_obj$f)
  tibble::tibble(
    Hypothesis      = label,
    Model           = sprintf("One-way: %s", group),
    k               = k,
    n_h             = round(n_h, 2),
    eta2            = round(eta2, 4),
    f_obs           = round(f_obs, 3),
    power_pct       = round(100 * pw, 1),
    f_MDES_80       = round(f_mdes, 3),
    d_approx_MDES   = round(2 * f_mdes, 3),
    df_num          = NA_real_,
    df_den          = NA_real_,
    eta2_partial    = NA_real_
  )
}

power_interaction <- function(a3, focal_term, cell_data, cell_factors, label) {
  rn <- rownames(a3)
  row_id <- match(focal_term, rn)
  if (is.na(row_id)) stop(sprintf("Focal term '%s' not found.", focal_term))
  Fval <- as.numeric(a3[row_id, "F value"])
  u    <- as.numeric(a3[row_id, "Df"])
  v    <- as.numeric(a3[nrow(a3), "Df"])
  eta2p <- (Fval * u) / (Fval * u + v)
  f_obs <- sqrt(eta2p / (1 - eta2p))
  f2    <- f_obs^2
  pw       <- pwr.f2.test(u = u, v = v, f2 = f2, sig.level = alpha)$power
  mdes_obj <- pwr.f2.test(u = u, v = v, sig.level = alpha, power = 0.80)
  f_mdes   <- sqrt(as.numeric(mdes_obj$f2))
  n_h <- NA_real_
  if (!is.null(cell_data) && !is.null(cell_factors)) {
    dcell <- cell_data %>% dplyr::mutate(.cell = interaction(!!!rlang::syms(cell_factors), drop = TRUE))
    ns    <- table(dcell$.cell)
    n_h   <- round(harmonic_mean(as.numeric(ns)), 2)
  }
  tibble::tibble(
    Hypothesis      = label,
    Model           = sprintf("LM focal: %s", focal_term),
    k               = NA_integer_,
    n_h             = n_h,
    eta2            = NA_real_,
    f_obs           = round(f_obs, 3),
    power_pct       = round(100 * pw, 1),
    f_MDES_80       = round(f_mdes, 3),
    d_approx_MDES   = round(2 * f_mdes, 3),
    df_num          = u,
    df_den          = v,
    eta2_partial    = round(eta2p, 4)
  )
}

res_H1 <- power_oneway(hypotheses_partisans, "response", "party_condition", "H1: party_condition main effect")

d_H2_pw <- hypotheses_partisans %>%
  dplyr::filter(party_condition %in% c("Inparty","Outparty")) %>%
  dplyr::mutate(party_condition = droplevels(party_condition))
res_H2 <- power_oneway(d_H2_pw, "response", "party_condition", "H2: Outparty vs Inparty")

res_H3 <- power_oneway(hypotheses_partisans, "response", "disruption_severity", "H3: severity main effect")

d_H5 <- hypotheses_partisans %>%
  filter(party_condition %in% c("Inparty","Outparty")) %>%
  mutate(
    party_condition = factor(party_condition, levels = c("Inparty","Outparty")),
    disruption_severity = factor(disruption_severity, levels = c("Mild","Severe"))
  )

# Fit model for coefficients (uses treatment contrasts by default: Inparty = reference)
m_H5 <- lm(response ~ party_condition * disruption_severity, data = d_H5)
cat("\nH5 OLS summary:\n"); print(summary(m_H5))

# Type III ANOVA with sum contrasts
old_contr_H5 <- options(contrasts = c("contr.sum","contr.poly"))
m_H5_anova <- lm(response ~ party_condition * disruption_severity, data = d_H5)
a3_H5 <- car::Anova(m_H5_anova, type = "III")
options(old_contr_H5)
cat("\nH5 Type-III ANOVA:\n"); print(a3_H5)

res_H5 <- power_interaction(
  a3 = a3_H5,
  focal_term = "party_condition:disruption_severity",
  cell_data = d_H5,
  cell_factors = c("party_condition","disruption_severity"),
  label = "H5: (In vs Out) x Severity interaction"
)

# NEW (safety): define d_H6 for H6 model (party_family subset)
d_H6 <- hypotheses_partisans %>%
  dplyr::filter(party_family %in% c("Democrat","Republican")) %>%
  dplyr::mutate(
    party_family = factor(party_family, levels = c("Democrat","Republican")),
    disruption_severity = factor(disruption_severity, levels = c("Mild","Severe"))
  )

# Fit model for coefficients (uses treatment contrasts by default: Democrat = reference)
m_H6 <- lm(response ~ party_family * disruption_severity, data = d_H6)
cat("\nH6 OLS summary:\n"); print(summary(m_H6))

# Type III ANOVA with sum contrasts
old_contr_H6 <- options(contrasts = c("contr.sum","contr.poly"))
m_H6_anova <- lm(response ~ party_family * disruption_severity, data = d_H6)
a3_H6 <- car::Anova(m_H6_anova, type = "III")
options(old_contr_H6)
cat("\nH6 Type-III ANOVA:\n"); print(a3_H6)

res_H6 <- power_interaction(
  a3 = a3_H6,
  focal_term = "party_family:disruption_severity",
  cell_data = d_H6,
  cell_factors = c("party_family","disruption_severity"),
  label = "H6: (Dem vs Rep) x Severity interaction"
)

# ============================================================
# emmeans simple-effects tables (A6b, A9b)
#   - A6b: H4 (formerly H5) Outparty–Inparty gaps within Mild/Severe
#   - A9b: H5 (formerly H6) Severe–Mild slopes within Democrat/Republican
# ============================================================

# ---- H4 (Inparty vs Outparty x Severity) -> Outparty–Inparty by severity
emm_H4 <- emmeans(m_H5, ~ party_condition * disruption_severity)
con_H4 <- pairs(emm_H4, by = "disruption_severity") %>% as.data.frame()

SI_Table_6b <- con_H4 %>%
  dplyr::mutate(
    `Disruption severity` = disruption_severity,
    Contrast = "Outparty – Inparty",
    Estimate = dplyr::if_else(grepl("Inparty - Outparty", contrast), -estimate, estimate),
    SE = SE,
    df = df,
    `t value` = abs(t.ratio),
    `p-value` = p.value
  ) %>%
  dplyr::select(`Disruption severity`, Contrast, Estimate, SE, df, `t value`, `p-value`) %>%
  dplyr::mutate(
    dplyr::across(c(Estimate, SE, `t value`), ~ round(., 3)),
    `p-value` = ifelse(`p-value` < 0.0001, "<0.0001", sprintf("%.4f", `p-value`))
  )

cat("\nA6b (H4) emmeans simple-effects table prepared.\n")
print(SI_Table_6b)

# ---- H5 (Democrat vs Republican x Severity) -> Severe–Mild within party
emm_H5 <- emmeans(m_H6, ~ party_family * disruption_severity)
con_H5 <- pairs(emm_H5, by = "party_family") %>% as.data.frame()

SI_Table_8b <- con_H5 %>%
  dplyr::mutate(
    Party = party_family,
    Contrast = "Severe – Mild",
    Estimate = dplyr::if_else(grepl("Mild - Severe", contrast), -estimate, estimate),
    SE = SE,
    df = df,
    `t value` = abs(t.ratio),
    `p-value` = p.value
  ) %>%
  dplyr::select(Party, Contrast, Estimate, SE, df, `t value`, `p-value`) %>%
  dplyr::mutate(
    dplyr::across(c(Estimate, SE, `t value`), ~ round(., 3)),
    `p-value` = ifelse(`p-value` < 0.0001, "<0.0001", sprintf("%.4f", `p-value`))
  )

cat("\nA9b (H5) emmeans simple-effects table prepared.\n")
print(SI_Table_8b)

power_summary <- dplyr::bind_rows(res_H1, res_H2, res_H3, res_H5, res_H6)
cat("\nPower summary (observed f, achieved power %, MDES @80%):\n")
print(power_summary)

cat("\nReporting sentences:\n")
report_sentence <- function(row) {
  sprintf(
    "Using the observed effect size (f = %.3f) and the harmonic mean group size (n_h = %.2f), post hoc power was estimated at %.1f%% (alpha = .05). A sensitivity analysis indicated that with the available sample size, the design had 80%% power to detect effects of at least f = %.3f (~ d - %.3f).",
    as.numeric(row$f_obs), as.numeric(row$n_h), as.numeric(row$power_pct),
    as.numeric(row$f_MDES_80), as.numeric(row$d_approx_MDES)
  )
}
for (i in seq_len(nrow(power_summary))) {
  r <- power_summary[i, ]
  cat(sprintf("[%s] %s\n", r$Hypothesis, report_sentence(r)))
}

d_rep <- hypotheses_partisans %>% filter(party_family == "Republican")
d_dem <- hypotheses_partisans %>% filter(party_family == "Democrat")

res_H1_R <- power_oneway(d_rep, "response", "party_condition", "H1 (Rep): party_condition main effect")
res_H1_D <- power_oneway(d_dem, "response", "party_condition", "H1 (Dem): party_condition main effect")
res_H2_R <- d_rep %>% filter(party_condition %in% c("Inparty","Outparty")) %>% mutate(party_condition = droplevels(party_condition)) %>% { power_oneway(., "response", "party_condition", "H2 (Rep): Outparty vs Inparty") }
res_H2_D <- d_dem %>% filter(party_condition %in% c("Inparty","Outparty")) %>% mutate(party_condition = droplevels(party_condition)) %>% { power_oneway(., "response", "party_condition", "H2 (Dem): Outparty vs Inparty") }
res_H3_R <- power_oneway(d_rep, "response", "disruption_severity", "H3 (Rep): severity main effect")
res_H3_D <- power_oneway(d_dem, "response", "disruption_severity", "H3 (Dem): severity main effect")

power_summary_party <- bind_rows(res_H1_R, res_H1_D, res_H2_R, res_H2_D, res_H3_R, res_H3_D)
cat("\nSubgroup power summary (Republicans/Democrats for H1-H3):\n")
print(power_summary_party)

# ============================================================
# STEP VI. Appendix: Dark personality (merge cleaning + models)
#   - Harmonize variables for covariate-adjusted DARK models
#   - Build DARK1-DARK5 tables
# ============================================================

cat("\n--- Appendix: Dark personality (H4, H7) - merged cleaning + models ---\n")

# Harmonize variables for covariate-adjusted analyses
dp_fix <- hypotheses_partisans %>%
  transmute(
    response,
    exploitation = C1_Dark_personality,                                 # 1-7
    party_condition = factor(party_condition, levels = c("No Party","Inparty","Outparty")),
    disruption_severity = factor(disruption_severity, levels = c("Mild","Severe")),
    party_family = factor(party_family, levels = c("Democrat","Republican")),
    partisanship_raw = factor(partisanship, levels = c("Moderate","Strong")),
    age = suppressWarnings(as.numeric(age)),
    gender_raw = if (haven::is.labelled(Q_Gender)) haven::as_factor(Q_Gender, levels = "labels") else Q_Gender,
    income_raw = if (haven::is.labelled(profile_gross_household)) haven::as_factor(profile_gross_household, levels = "labels") else profile_gross_household,
    race_raw = if (haven::is.labelled(race)) haven::as_factor(race, levels = "labels") else race
  ) %>%
  filter(!is.na(response), !is.na(exploitation)) %>%
  mutate(
    gender_raw = as.character(gender_raw),
    income_raw = as.character(income_raw)
  ) %>%
  filter(
    gender_raw %in% c("Male","Female"),
    party_family %in% c("Democrat","Republican"),
    !income_raw %in% c("Don't know","Prefer not to answer")
  ) %>%
  mutate(
    gender = factor(gender_raw, levels = c("Male","Female")),              # ref = Male
    race = factor(ifelse(race_raw %in% c("White"), "White", "Non-White"),
                  levels = c("White","Non-White")),                        # ref = White
    partisanship = factor(ifelse(partisanship_raw == "Strong", "Strong", "Weak"),
                          levels = c("Weak","Strong")),                     # ref = Weak
    income_group = factor(
      dplyr::case_when(
        income_raw %in% c("Less than $10,000",
                          "$10,000 - $19,999",
                          "$20,000 - $29,999") ~ "Low",
        income_raw %in% c("$30,000 - $39,999", "$40,000 - $49,999",
                          "$50,000 - $59,999", "$60,000 - $69,999",
                          "$70,000 - $79,999", "$80,000 - $99,999") ~ "Middle",
        TRUE ~ "High"
      ),
      levels = c("Low","Middle","High")                                     # ref = Low
    ),
    # Convenience dummies
    female  = as.integer(gender == "Female"),
    nonwhite = as.integer(race == "Non-White"),
    party_out_vs_control = as.integer(party_condition == "Outparty"),
    party_in_vs_control  = as.integer(party_condition == "Inparty"),
    partisan_strong = as.integer(partisanship == "Strong"),
    disruption_severity_mild = as.integer(dplyr::recode(disruption_severity, Mild = 1L, Severe = 0L)),
    income_middle = as.integer(income_group == "Middle"),
    income_high   = as.integer(income_group == "High"),
    republican = as.integer(party_family == "Republican"),
    age10 = age/10
  )

cat(sprintf("Harmonized N (dp_fix) = %d\n", nrow(dp_fix)))

# --- Lock treatment baselines inside modeling data ---
dp_fix <- within(dp_fix, {
  party_condition       <- factor(party_condition,       levels = c("No Party","Inparty","Outparty"))
  disruption_severity   <- factor(disruption_severity,   levels = c("Mild","Severe"))
  party_family          <- factor(party_family,          levels = c("Democrat","Republican"))
  partisanship          <- factor(partisanship,          levels = c("Weak","Strong"))
  gender                <- factor(gender,                levels = c("Male","Female"))
  income_group          <- factor(income_group,          levels = c("Low","Middle","High"))
  race                  <- factor(race,                  levels = c("White","Non-White"))
})
contrasts(dp_fix$party_condition)     <- contr.treatment(nlevels(dp_fix$party_condition),     base = 1)
contrasts(dp_fix$disruption_severity) <- contr.treatment(nlevels(dp_fix$disruption_severity), base = 1)
contrasts(dp_fix$party_family)        <- contr.treatment(nlevels(dp_fix$party_family),        base = 1)
contrasts(dp_fix$partisanship)        <- contr.treatment(nlevels(dp_fix$partisanship),        base = 1)
contrasts(dp_fix$gender)              <- contr.treatment(nlevels(dp_fix$gender),              base = 1)
contrasts(dp_fix$income_group)        <- contr.treatment(nlevels(dp_fix$income_group),        base = 1)
contrasts(dp_fix$race)                <- contr.treatment(nlevels(dp_fix$race),                base = 1)

# H4 baseline (bivariate)
m_H4_biv <- lm(response ~ exploitation, data = dp_fix)

# H4 covariate-adjusted (continuous exploitation)
m_H4_ols <- lm(
  response ~ exploitation + party_condition + disruption_severity +
    partisanship + party_family + age10 + gender + income_group + race,
  data = dp_fix
)

# H4 alternative: High exploitation (6-7)
dp_fix <- dp_fix %>% mutate(exploit_hi_67 = as.integer(exploitation >= 6))
m_H4_ols67 <- lm(
  response ~ exploit_hi_67 + party_condition + disruption_severity +
    partisanship + party_family + age10 + gender + income_group + race,
  data = dp_fix
)

# H4 standardized DV & exploitation (report std beta)
m_H4_ols_std <- lm(
  scale(response) ~ scale(exploitation) + party_condition + disruption_severity +
    partisanship + party_family + scale(age10) + gender + income_group + race,
  data = dp_fix
)

# H7: GAM (smooth of exploitation) with covariates
gam_H7 <- mgcv::gam(
  response ~ s(exploitation, k = 4) + party_condition + disruption_severity +
    partisanship + party_family + s(age10) + gender + income_group + race,
  data = dp_fix
)

# Tidy DARK tables
tidy_keep <- function(m){
  broom::tidy(m, conf.int = TRUE) %>%
    dplyr::mutate(across(where(is.numeric), ~round(., 3))) %>%
    dplyr::rename(`CI lower` = conf.low, `CI upper` = conf.high,
                  `Std. Error` = std.error, `t value` = statistic, `p-value` = p.value)
}
DARK1 <- tidy_keep(m_H4_biv)
DARK2 <- tidy_keep(m_H4_ols)
DARK3 <- tidy_keep(m_H4_ols67)
DARK4 <- tidy_keep(m_H4_ols_std)
DARK5 <- broom::tidy(gam_H7) %>% dplyr::mutate(across(where(is.numeric), ~round(., 3)))

cat("\n(H4/H7) DARK tables prepared: DARK1-DARK5\n")

# ============================================================
# STEP VII. Robustness with covariates (H1-H6) - build R0-R6
# ============================================================

cat("\n--- Robustness: H1-H6 OLS with covariates (R#) ---\n")

# R0: Baseline OLS (party only) using dp_fix (treatment baseline = No Party)
R0_base <- lm(response ~ party_condition, data = dp_fix)
R0_tbl_coef <- tidy_keep(R0_base)

# H1: party_condition main effect (baseline = No Party) + covariates
H1_ols_cov <- lm(response ~ party_condition + partisanship + disruption_severity +
                   party_family + age10 + gender + income_group + race,
                 data = dp_fix)
cat("\n--- H1 OLS with covariates (baseline = No Party) ---\n"); print(summary(H1_ols_cov))

# H2: Outparty vs Inparty (subset; baseline = Inparty)
d_H2_cov <- dp_fix %>%
  dplyr::filter(party_condition %in% c("Inparty","Outparty")) %>%
  dplyr::mutate(party_condition = droplevels(party_condition),
                party_condition = relevel(party_condition, ref = "Inparty"))
contrasts(d_H2_cov$party_condition) <- contr.treatment(2, base = 1)

H2_ols_cov <- lm(response ~ party_condition + partisanship + disruption_severity +
                   party_family + age10 + gender + income_group + race,
                 data = d_H2_cov)

# H3: Severity main effect (baseline = Mild)
H3_ols_cov <- lm(response ~ disruption_severity + partisanship + party_condition +
                   party_family + age10 + gender + income_group + race,
                 data = dp_fix)

# H5: (In vs Out) - Severity
d_H5_cov <- dp_fix %>%
  dplyr::filter(party_condition %in% c("Inparty","Outparty")) %>%
  dplyr::mutate(
    party_condition = droplevels(party_condition),
    party_condition = relevel(party_condition, ref = "Inparty"),
    disruption_severity = factor(disruption_severity, levels = c("Mild","Severe"))
  )
contrasts(d_H5_cov$party_condition)     <- contr.treatment(2, base = 1)
contrasts(d_H5_cov$disruption_severity) <- contr.treatment(2, base = 1)

H5_ols_cov <- lm(response ~ party_condition * disruption_severity + partisanship +
                   party_family + age10 + gender + income_group + race,
                 data = d_H5_cov)

# H6: party_family - Severity (baseline = Democrat, Mild)
d_H6_cov <- dp_fix %>%
  dplyr::filter(party_family %in% c("Democrat","Republican")) %>%
  dplyr::mutate(disruption_severity = factor(disruption_severity, levels = c("Mild","Severe")))
contrasts(d_H6_cov$party_family)        <- contr.treatment(2, base = 1)
contrasts(d_H6_cov$disruption_severity) <- contr.treatment(2, base = 1)

H6_ols_cov <- lm(response ~ party_family * disruption_severity + partisanship +
                   party_condition + age10 + gender + income_group + race,
                 data = d_H6_cov)

# Type III ANOVAs for H5 & H6 (with sum contrasts only for ANOVA table)
old_contr_rob <- options(contrasts = c("contr.sum","contr.poly"))
a3_H5_cov <- car::Anova(H5_ols_cov, type = "III")
a3_H6_cov <- car::Anova(H6_ols_cov, type = "III")
options(old_contr_rob)

# Tidy R# tables
R1_H1 <- tidy_keep(H1_ols_cov)
R2_H2 <- tidy_keep(H2_ols_cov)
R3_H3 <- tidy_keep(H3_ols_cov)
R4_H5 <- tidy_keep(H5_ols_cov)
R5_H6 <- tidy_keep(H6_ols_cov)

R6_fit <- dplyr::bind_rows(
  dplyr::mutate(broom::glance(R0_base),     Model = "R0: party only (No / In / Out)"),
  dplyr::mutate(broom::glance(H1_ols_cov), Model = "H1: party_condition + covariates"),
  dplyr::mutate(broom::glance(H2_ols_cov), Model = "H2: Outparty vs Inparty + covariates"),
  dplyr::mutate(broom::glance(H3_ols_cov), Model = "H3: severity + covariates"),
  dplyr::mutate(broom::glance(H5_ols_cov), Model = "H5: (In/Out) - severity + covariates"),
  dplyr::mutate(broom::glance(H6_ols_cov), Model = "H6: party_family - severity + covariates")
) %>%
  dplyr::select(Model, nobs, r.squared, adj.r.squared, sigma, statistic, df, df.residual, p.value) %>%
  dplyr::mutate(across(where(is.numeric), ~round(., 3))) %>%
  dplyr::rename(`N`=nobs, `R^2`=r.squared, `Adj. R^2`=adj.r.squared, `Residual SE`=sigma,
                `Model F`=statistic, `Model df1`=df, `Model df2`=df.residual, `Model p`=p.value)

cat("\nRobustness R tables prepared: R0-R6\n")






###
# ============================================================
# STEP Y. Build all tables for export (D, A, E, P)
#   - Creates Descr_Table_* , SI_Table_* , E- and P-tables
#   - Uses objects defined in earlier steps; re-computes if needed
# ============================================================

# ---- Descriptive tables (D#) ----
# D1
Descr_Table_1 <- tibble::tibble(
  Metric = c("Total with non-NA response",
             "Inattentive (QAttention != 4)",
             "Non-partisans (incl. true independents)",
             "Retained (attentive partisans)"),
  N = c(N_all, N_inattentive, N_nonpartisans, N_retained)
)

# D2
Descr_Table_2 <- tab_npct(hypotheses_partisans$partisanship)

# D3
Descr_Table_3 <- tab_npct(hypotheses_partisans$party_family)

# D4
gender_vecY <- if (haven::is.labelled(hypotheses_partisans$Q_Gender)) {
  haven::as_factor(hypotheses_partisans$Q_Gender, levels = "labels")
} else hypotheses_partisans$Q_Gender
Descr_Table_4 <- tab_npct(gender_vecY)

# D5A / D5B
age_sum <- summary(hypotheses_partisans$age)
Descr_Table_5A <- tibble::tibble(
  Statistic = names(age_sum),
  Value = as.numeric(age_sum)
)
age_binsY <- cut(hypotheses_partisans$age,
                 breaks = c(18,25,35,45,55,65,75,85,Inf), right = FALSE)
Descr_Table_5B <- tab_npct(age_binsY)

# D6
income_cleanY <- if (haven::is.labelled(hypotheses_partisans$profile_gross_household)) {
  haven::as_factor(hypotheses_partisans$profile_gross_household, levels = "labels")
} else hypotheses_partisans$profile_gross_household
if (is.factor(income_cleanY)) {
  income_cleanY <- droplevels(income_cleanY[income_cleanY != "Don't know/Prefer not to say"])
} else {
  income_cleanY[income_cleanY == 99] <- NA
}
Descr_Table_6 <- tab_npct(income_cleanY)

# D7
race_vecY <- if (haven::is.labelled(hypotheses_partisans$race)) {
  haven::as_factor(hypotheses_partisans$race, levels = "labels")
} else hypotheses_partisans$race
Descr_Table_7 <- tab_npct(race_vecY)

# D8–D10
Descr_Table_8  <- hypotheses_partisans %>% count(party_condition, name = "N")
Descr_Table_9  <- hypotheses_partisans %>% count(disruption_severity, name = "N")
Descr_Table_10 <- hypotheses_partisans %>%
  count(party_condition, disruption_severity, name = "N") %>%
  arrange(party_condition, disruption_severity)

# ---- Main A-tables (SI_Table_*) ----
# A1: Means & 95% CIs by group (Overall, Republicans, Democrats) grouped by party_condition × disruption_severity 

mk_sum_tbl <- function(df, group_label) {
  df %>%
    dplyr::group_by(party_condition, disruption_severity) %>%
    dplyr::summarise(
      N  = dplyr::n(),
      Mean = mean(response, na.rm = TRUE),
      SD   = sd(response, na.rm = TRUE),
      SE   = SD / sqrt(N),
      CI_lower = Mean - qt(0.975, df = N - 1) * SE,
      CI_upper = Mean + qt(0.975, df = N - 1) * SE,
      .groups = "drop"
    ) %>%
    dplyr::mutate(Group = group_label)
}

SI_Table_2 <- dplyr::bind_rows(
  mk_sum_tbl(hypotheses_partisans, "All"),
  mk_sum_tbl(dplyr::filter(hypotheses_partisans, party_family == "Republican"), "Republicans"),
  mk_sum_tbl(dplyr::filter(hypotheses_partisans, party_family == "Democrat"),   "Democrats")
) %>%
  dplyr::mutate(
    dplyr::across(c(Mean, SD, CI_lower, CI_upper), ~ round(., 2)),
    `95% CI [lower–upper]` = sprintf("[%.2f, %.2f]", CI_lower, CI_upper)
  ) %>%
  dplyr::rename(
    `Party condition` = party_condition,
    `Disruption severity` = disruption_severity
  ) %>%
  dplyr::select(Group, `Party condition`, `Disruption severity`, N, Mean, SD, `95% CI [lower–upper]`)


# Ensure Welch tests/ANOVAs exist (recompute if missing)
if (!exists("t_H1_out_vs_no")) {
  no_party  <- filter(hypotheses_partisans, party_condition == "No Party")
  inparty   <- filter(hypotheses_partisans, party_condition == "Inparty")
  outparty  <- filter(hypotheses_partisans, party_condition == "Outparty")
  t_H1_out_vs_no <- t.test(outparty$response, no_party$response, var.equal = FALSE)
  t_H1_in_vs_no  <- t.test(inparty$response,  no_party$response, var.equal = FALSE)
  t_H2           <- t.test(outparty$response, inparty$response,  alternative="greater", var.equal = FALSE)
}
# A2: H1 (two-sided Welch) - WITH CI BOUNDS AND GROUP MEANS
mk_welch_row <- function(ttest_obj, comp_label, mean1, mean2) {
  tibble::tibble(
    Comparison = comp_label,
    `T-value` = unname(ttest_obj$statistic),
    df = unname(ttest_obj$parameter),
    `P-value` = unname(ttest_obj$p.value),
    `CI,lower` = ttest_obj$conf.int[1],
    `CI,upper` = ttest_obj$conf.int[2],
    `Mean(NoParty)` = mean1,
    `Mean(Outparty/Inparty)` = mean2
  )
}

SI_Table_3 <- bind_rows(
  mk_welch_row(
    t_H1_out_vs_no, 
    "No Party vs Outparty",
    mean(no_party$response, na.rm=TRUE),
    mean(outparty$response, na.rm=TRUE)
  ),
  mk_welch_row(
    t_H1_in_vs_no,
    "No Party vs Inparty", 
    mean(no_party$response, na.rm=TRUE),
    mean(inparty$response, na.rm=TRUE)
  )
) %>% 
  mutate(across(where(is.numeric), ~round(., 3)))

# A3: H2 (one-sided Welch Outparty > Inparty) - WITH CI BOUNDS AND GROUP MEANS
SI_Table_4 <- tibble::tibble(
  Comparison = "Outparty vs Inparty (one-sided)",
  `T-value` = unname(t_H2$statistic),
  df = unname(t_H2$parameter),
  `P-value` = unname(t_H2$p.value),
  `CI,lower` = t_H2$conf.int[1],
  `CI,upper` = t_H2$conf.int[2],
  `Mean(Outparty)` = mean(outparty$response, na.rm=TRUE),
  `Mean(Inparty)` = mean(inparty$response, na.rm=TRUE)
) %>% 
  mutate(across(where(is.numeric), ~round(., 3)))

# A4: H3 (overall & within party_condition) one-sided (Severe > Mild) - WITH CI BOUNDS AND GROUP MEANS
mk_t_row <- function(ttest_obj, sample_label, mean_mild, mean_severe) {
  tibble::tibble(
    Sample = sample_label,
    `T-value` = unname(ttest_obj$statistic),
    df = unname(ttest_obj$parameter),
    `P-value` = unname(ttest_obj$p.value),
    `CI,lower` = ttest_obj$conf.int[1],
    `CI,upper` = ttest_obj$conf.int[2],
    `Mean(Mild)` = mean_mild,
    `Mean(Severe)` = mean_severe
  )
}

mild_all   <- filter(hypotheses_partisans, disruption_severity=="Mild")$response
severe_all <- filter(hypotheses_partisans, disruption_severity=="Severe")$response
mild_no    <- filter(no_party,  disruption_severity=="Mild")$response
severe_no  <- filter(no_party,  disruption_severity=="Severe")$response
mild_out   <- filter(outparty,  disruption_severity=="Mild")$response
severe_out <- filter(outparty,  disruption_severity=="Severe")$response
mild_in    <- filter(inparty,   disruption_severity=="Mild")$response
severe_in  <- filter(inparty,   disruption_severity=="Severe")$response

SI_Table_5 <- bind_rows(
  mk_t_row(t_H3_all, "All", 
           mean(mild_all, na.rm=TRUE), 
           mean(severe_all, na.rm=TRUE)),
  mk_t_row(t_H3_no, "No Party",
           mean(mild_no, na.rm=TRUE),
           mean(severe_no, na.rm=TRUE)),
  mk_t_row(t_H3_out, "Outparty",
           mean(mild_out, na.rm=TRUE),
           mean(severe_out, na.rm=TRUE)),
  mk_t_row(t_H3_in, "Inparty",
           mean(mild_in, na.rm=TRUE),
           mean(severe_in, na.rm=TRUE))
) %>% 
  mutate(across(where(is.numeric), ~round(., 3)))

# A5: Omnibus ANOVAs & 2x3 factorial (Type III)
if (!exists("a_H1")) a_H1 <- aov(response ~ party_condition, data = hypotheses_partisans)
if (!exists("a_H2")) a_H2 <- aov(response ~ factor(party_condition, levels=c("Inparty","Outparty")),
                                 data = filter(hypotheses_partisans, party_condition %in% c("Inparty","Outparty")))
if (!exists("a_H3")) a_H3 <- aov(response ~ disruption_severity, data = hypotheses_partisans)
if (!exists("a3_2x3")) {
  op_contrY <- options(contrasts = c("contr.sum","contr.poly"))
  m_2x3Y <- lm(response ~ party_condition * disruption_severity, data = hypotheses_partisans)
  a3_2x3 <- car::Anova(m_2x3Y, type = "III")
  options(op_contrY)
}
grab_aov <- function(fit, term_name) {
  s <- summary(fit)[[1]]
  tibble::tibble(
    Term = term_name,
    df1 = as.numeric(s[1,"Df"]),
    df2 = as.numeric(s[nrow(s),"Df"]),
    F = as.numeric(s[1,"F value"]),
    p_value = as.numeric(s[1,"Pr(>F)"])
  )
}
grab_typeIII <- function(a3) {
  as.data.frame(a3) %>%
    tibble::rownames_to_column("Term") %>%
    dplyr::rename(df1 = Df, F = `F value`, p_value = `Pr(>F)`) %>%
    mutate(df2 = a3[nrow(a3), "Df"]) %>%
    select(Term, df1, df2, F, p_value)
}
SI_Table_2A <- bind_rows(
  grab_aov(a_H1, "party_condition (H1)"),
  grab_aov(a_H2, "Inparty vs Outparty (H2)"),
  grab_aov(a_H3, "disruption_severity (H3)"),
  grab_typeIII(a3_2x3) %>% mutate(Term = paste0("Type III: ", Term))
) %>% mutate(across(where(is.numeric), ~round(., 4)))

# A6–A11 from models already fit (m_H5 / a3_H5 / m_H6 / a3_H6)
SI_Table_6  <- broom::tidy(m_H5, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~round(., 3))) %>%
  rename(`Std. Error`=std.error, `t value`=statistic, `p-value`=p.value,
         `CI lower`=conf.low, `CI upper`=conf.high)

SI_Table_6B <- broom::glance(m_H5) %>%
  mutate(across(where(is.numeric), ~round(., 3)))

SI_Table_7  <- as.data.frame(a3_H5) %>%
  tibble::rownames_to_column("Term") %>%
  rename(df = Df, `F value`=`F value`, `p-value`=`Pr(>F)`) %>%
  mutate(across(where(is.numeric), ~round(., 4)))

SI_Table_8  <- broom::tidy(m_H6, conf.int = TRUE) %>%
  mutate(across(where(is.numeric), ~round(., 3))) %>%
  rename(`Std. Error`=std.error, `t value`=statistic, `p-value`=p.value,
         `CI lower`=conf.low, `CI upper`=conf.high)

SI_Table_8B <- broom::glance(m_H6) %>%
  mutate(across(where(is.numeric), ~round(., 3)))

SI_Table_9  <- as.data.frame(a3_H6) %>%
  tibble::rownames_to_column("Term") %>%
  rename(df = Df, `F value`=`F value`, `p-value`=`Pr(>F)`) %>%
  mutate(across(where(is.numeric), ~round(., 4)))

# (A6b and A9b were already created above as SI_Table_6b and SI_Table_8b)

# ---- Exploratory chi-square tables (E#) ----
# Reuse G and prop_diff_stats defined earlier
build_H2_chi <- function() {
  rows <- lapply(names(G), function(lbl){
    dd <- hypotheses_partisans[G[[lbl]](hypotheses_partisans), , drop = FALSE]
    st <- prop_diff_stats(dd, dd$party_condition=="Outparty", dd$party_condition=="Inparty", alt="greater")
    tibble::tibble(Group = lbl, est = st$est, p_value = st$p, CI_low = st$ci[1], CI_high = st$ci[2],
                   n_Out = st$n[1], n_In = st$n[2])
  })
  bind_rows(rows) %>%
    mutate(across(where(is.numeric), ~round(., 6)))
}
SI_Table_10_H2_Chi <- build_H2_chi()

build_H1_chi <- function() {
  rows <- lapply(names(G), function(lbl){
    dd <- hypotheses_partisans[G[[lbl]](hypotheses_partisans), , drop = FALSE]
    st1 <- prop_diff_stats(dd, dd$party_condition=="No Party", dd$party_condition=="Outparty", alt="less")
    st2 <- prop_diff_stats(dd, dd$party_condition=="No Party", dd$party_condition=="Inparty",  alt="greater")
    tibble::tibble(
      Group = lbl,
      Comparison = c("NoParty - Outparty (one-sided, less)", "NoParty - Inparty (one-sided, greater)"),
      est = c(st1$est, st2$est),
      p_value = c(st1$p, st2$p),
      CI_low = c(st1$ci[1], st2$ci[1]),
      CI_high = c(st1$ci[2], st2$ci[2]),
      n_No = sum(dd$party_condition=="No Party"),
      n_Out = sum(dd$party_condition=="Outparty"),
      n_In  = sum(dd$party_condition=="Inparty")
    )
  })
  bind_rows(rows) %>% mutate(across(where(is.numeric), ~round(., 6)))
}
SI_Table_11_H1_Chi <- build_H1_chi()

build_H3_chi <- function() {
  rows <- lapply(names(G), function(lbl){
    dd <- hypotheses_partisans[G[[lbl]](hypotheses_partisans), , drop = FALSE]
    st <- prop_diff_stats(dd, dd$disruption_severity=="Severe", dd$disruption_severity=="Mild", alt="greater")
    tibble::tibble(Group = lbl, est = st$est, p_value = st$p, CI_low = st$ci[1], CI_high = st$ci[2],
                   n_Severe = st$n[1], n_Mild = st$n[2])
  })
  bind_rows(rows) %>%
    mutate(across(where(is.numeric), ~round(., 6)))
}
SI_Table_12_H3_Chi <- build_H3_chi()

# ---- Post hoc power/sensitivity tables (P#) ----
power_tbl_export <- power_summary

SI_Table_10B <- tibble::tibble(
  Hypothesis = power_summary$Hypothesis,
  Sentence = vapply(seq_len(nrow(power_summary)), function(i) {
    r <- power_summary[i, ]
    sprintf(
      "Using the observed effect size (f = %.3f) and the harmonic mean group size (n_h = %.2f), post hoc power was estimated at %.1f%% (alpha = .05). A sensitivity analysis indicated that with the available sample size, the design had 80%% power to detect effects of at least f = %.3f (~ d - %.3f).",
      as.numeric(r$f_obs), as.numeric(r$n_h), as.numeric(r$power_pct),
      as.numeric(r$f_MDES_80), as.numeric(r$d_approx_MDES)
    )
  }, character(1))
)

tbl_party_export <- power_summary_party




# ============================================================
# STEP Z. Export ONE Word file with all tables:
#   Descriptives (D#), Main analyses (A#), Exploratory (E#),
#   Post hoc power (P#), Robustness (R#), Dark personality (DARK#)
# ============================================================

need <- c("dplyr","tibble","broom","officer","flextable")
to_install <- need[!sapply(need, requireNamespace, quietly = TRUE)]
if (length(to_install)) install.packages(to_install)
library(dplyr); library(tibble); library(broom); library(officer); library(flextable)

to_safe_utf8 <- function(s) {
  s <- as.character(s)
  s <- iconv(s, from = "", to = "UTF-8", sub = "")
  s <- gsub("\u00D7","x", s, useBytes = FALSE)
  s <- gsub("[\u2012\u2013\u2014\u2015\u2212]","-", s, useBytes = FALSE)
  s <- gsub("\u00A0"," ", s, useBytes = FALSE)
  enc2utf8(s)
}
sanitize_df <- function(df) {
  nm <- names(df); nm <- to_safe_utf8(nm); names(df) <- nm
  for (j in seq_along(df)) {
    x <- df[[j]]
    if (is.factor(x)) {
      lev <- to_safe_utf8(levels(x))
      x <- factor(to_safe_utf8(as.character(x)), levels = to_safe_utf8(lev))
    } else if (is.character(x)) {
      x <- to_safe_utf8(x)
    }
    df[[j]] <- x
  }
  df
}
add_tbl <- function(doc, title, x, hstyle = NULL) {
  if (is.null(hstyle)) hstyle <- if ("heading 2" %in% styles_info(doc)$style_name) "heading 2" else "Heading 2"
  title <- to_safe_utf8(title)
  x <- sanitize_df(as.data.frame(x))
  doc <- body_add_par(doc, value = title, style = hstyle)
  ft  <- flextable::regulartable(x)
  ft  <- flextable::autofit(ft)
  flextable::body_add_flextable(doc, ft)
}

# ---- SAFE EXPORT HELPERS (NEW) ----
add_tbl_safe <- function(doc, title, obj_name, hstyle = NULL) {
  if (is.null(hstyle)) hstyle <- if ("heading 2" %in% styles_info(doc)$style_name) "heading 2" else "Heading 2"
  if (!exists(obj_name, inherits = TRUE)) {
    message(sprintf("Skipping '%s' (%s): object not found.", title, obj_name))
    return(doc)
  }
  x <- get(obj_name, inherits = TRUE)
  x <- sanitize_df(as.data.frame(x))
  title <- to_safe_utf8(title)
  doc <- body_add_par(doc, value = title, style = hstyle)
  ft  <- flextable::regulartable(x)
  ft  <- flextable::autofit(ft)
  flextable::body_add_flextable(doc, ft)
}
add_section <- function(doc, title, hstyle) {
  body_add_par(doc, value = title, style = hstyle)
}

# ---- LISTS OF TABLES TO EXPORT (name -> title) ----
D_tables <- list(
  Descr_Table_1  = "D1. Sample sizes and retention",
  Descr_Table_2  = "D2. Partisanship among retained (Strong/Moderate)",
  Descr_Table_3  = "D3. Party family among retained (Democrat/Republican)",
  Descr_Table_4  = "D4. Gender distribution among retained",
  Descr_Table_5A = "D5A. Age summary among retained",
  Descr_Table_5B = "D5B. Age distribution (10-year bins) among retained",
  Descr_Table_6  = "D6. Income distribution among retained",
  Descr_Table_7  = "D7. Race/ethnicity among retained",
  Descr_Table_8  = "D8. Cell counts by party condition",
  Descr_Table_9  = "D9. Cell counts by disruption severity",
  Descr_Table_10 = "D10. Cell counts by party condition - disruption severity"
)
A_tables <- list(
  SI_Table_2  = "A1. Means and 95% CIs by group (Overall, Republican, Democrat)",
  SI_Table_3  = "A2. H1 Welch: No Party vs Outparty / Inparty (two-sided)",
  SI_Table_4  = "A3. H2 Welch: Outparty > Inparty (one-sided)",
  SI_Table_5  = "A4. H3 Welch: Severe > Mild (overall & within party)",
  SI_Table_2A = "A5. Omnibus ANOVAs and 2-3 factorial (Type III)",
  SI_Table_6  = "A6. H5 OLS coefficients with 95% CIs (Inparty vs Outparty - Severity)",
  SI_Table_6b = "A6b. H4 simple effects (Outparty–Inparty by severity, emmeans)",  # NEW
  SI_Table_6B = "A7. H5 model fit statistics",
  SI_Table_7  = "A8. H5 Type-III ANOVA",
  SI_Table_8  = "A9. H6 OLS coefficients with 95% CIs (Democrat vs Republican - Severity)",
  SI_Table_8b = "A9b. H5 simple effects (Severe–Mild by party, emmeans)",          # NEW
  SI_Table_8B = "A10. H6 model fit statistics",
  SI_Table_9  = "A11. H6 Type-III ANOVA"
)
E_tables <- list(
  SI_Table_10_H2_Chi = "E1. H2: One-tailed Chi-square (Outparty - Inparty)",
  SI_Table_11_H1_Chi = "E2. H1: One-tailed Chi-square (No Party - Outparty; No Party - Inparty)",
  SI_Table_12_H3_Chi = "E3. H3: One-tailed Chi-square (Severe - Mild)"
)
P_tables <- list(
  power_tbl_export = "P1. Post hoc power and MDES summary (H1-H6)",
  SI_Table_10B     = "P2. Reporting sentences (power/MDES)",
  tbl_party_export = "P3. Post hoc power and MDES by party subgroup (H1-H3)"
)
R_tables <- list(
  R0_tbl_coef = "R0. H1 baseline OLS (party only; ref = No Party)",
  R1_H1       = "R1. H1 OLS + covariates",
  R2_H2       = "R2. H2 OLS + covariates (In vs Out)",
  R3_H3       = "R3. H3 OLS + covariates (Mild vs Severe)",
  R4_H5       = "R4. H5 OLS + covariates (In/Out - Severity)",
  R5_H6       = "R5. H6 OLS + covariates (Party family - Severity)",
  R6_fit      = "R6. Model fit summary (R0-H6)"
)
DARK_tables <- list(
  DARK1 = "DARK1. H4 bivariate OLS (Exploitation - Response)",
  DARK2 = "DARK2. H4 covariate-adjusted OLS (Exploitation, continuous)",
  DARK3 = "DARK3. H4 covariate-adjusted OLS (High exploitation 6-7)",
  DARK4 = "DARK4. H4 standardized OLS (Std. DV & Exploitation)",
  DARK5 = "DARK5. H7 GAM summary (smooth of exploitation, covariate-adjusted)"
)

# ---- Build the Word doc (SAFE EXPORT) ----
doc <- read_docx()
styles <- styles_info(doc)
h1 <- if ("heading 1" %in% styles$style_name) "heading 1" else "Heading 1"
h2 <- if ("heading 2" %in% styles$style_name) "heading 2" else "Heading 2"

# Title
doc <- body_add_par(doc, value = to_safe_utf8("Appendix: All Tables (D, A, E, P, R, DARK)"), style = h1)

# D tables
doc <- add_section(doc, "Descriptives (D-tables)", h1)
for (nm in names(D_tables)) doc <- add_tbl_safe(doc, D_tables[[nm]], nm, h2)

# A tables
doc <- add_section(doc, "Main analyses (A-tables)", h1)
for (nm in names(A_tables)) doc <- add_tbl_safe(doc, A_tables[[nm]], nm, h2)

# E tables
doc <- add_section(doc, "Exploratory analyses (E-tables)", h1)
for (nm in names(E_tables)) doc <- add_tbl_safe(doc, E_tables[[nm]], nm, h2)

# P tables
doc <- add_section(doc, "Post hoc power (P-tables)", h1)
for (nm in names(P_tables)) doc <- add_tbl_safe(doc, P_tables[[nm]], nm, h2)

# R tables
doc <- add_section(doc, "Robustness: OLS with covariates (R-tables)", h1)
for (nm in names(R_tables)) doc <- add_tbl_safe(doc, R_tables[[nm]], nm, h2)

# DARK tables
doc <- add_section(doc, "Dark personality appendix (DARK-tables)", h1)
for (nm in names(DARK_tables)) doc <- add_tbl_safe(doc, DARK_tables[[nm]], nm, h2)

# Save
outfile_all <- sprintf("Appendix_ALL_Tables_%s.docx", format(Sys.Date(), "%Y-%m-%d"))
print(doc, target = outfile_all)
cat("\nSaved Word file:", outfile_all, "\n")

# ============================================================
# STEP F. Figures for Main Paper (H1-H3) - attentive partisans only
# ============================================================

library(dplyr)
library(ggplot2)
library(patchwork)

pd <- position_dodge(width = 0.7)

# Figure 1 - H1
color_red  <- "#E34A33"  # Republicans
color_blue <- "#4F6DB8"  # Democrats
color_gray <- "gray60"   # Average

mk_sum <- function(df, lab){
  df %>%
    group_by(party_condition) %>%
    summarise(mean = mean(response, na.rm=TRUE),
              sd   = sd(response, na.rm=TRUE),
              n    = n(), .groups="drop") %>%
    mutate(group_type = lab,
           se = sd/sqrt(n),
           ci95_low = mean - 1.96*se,
           ci95_high= mean + 1.96*se,
           ci83_low = mean - 1.39*se,
           ci83_high= mean + 1.39*se)
}

d_rep <- mk_sum(hypotheses_partisans %>% filter(party_family=="Republican"), "Republicans")
d_dem <- mk_sum(hypotheses_partisans %>% filter(party_family=="Democrat"),   "Democrats")
d_avg <- mk_sum(hypotheses_partisans,                                        "Average")

d1_fig <- bind_rows(d_rep, d_avg, d_dem) %>%
  mutate(
    group_type = factor(group_type, levels=c("Republicans","Average","Democrats")),
    party_condition = factor(party_condition, levels=c("No Party","Inparty","Outparty"))
  )

Figure1 <- ggplot(d1_fig, aes(x = party_condition, y = mean, color = group_type)) +
  geom_point(position = pd, size = 4) +
  geom_errorbar(aes(ymin = ci83_low, ymax = ci83_high),
                position = pd, width = 0, size = 2, alpha = 0.8) +
  geom_errorbar(aes(ymin = ci95_low, ymax = ci95_high),
                position = pd, width = 0, size = 0.6, alpha = 0.8) +
  scale_color_manual(values = c("Republicans" = color_red,
                                "Average"     = color_gray,
                                "Democrats"   = color_blue)) +
  labs(x = "Victim Identity Marker",
       y = "Violence justification (full scale 1-7)",
       color = "Group") +
  scale_y_continuous(limits = c(2, 3.5)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.border     = element_rect(color = "black", fill = NA, linewidth = 0.5),
        axis.text.x      = element_text(size = 12, color = "black"),
        axis.text.y      = element_text(size = 12, color = "black"),
        axis.title       = element_text(color = "black"))
Figure1
ggsave("Figure_1_H1_means.png", Figure1, width=7.5, height=4.5, dpi=300)

# Figure 2 - H2: proportion gaps
hp <- hypotheses_partisans %>%
  mutate(
    condone = as.numeric(response) >= 5,
    part_strength = case_when(
      partisanship == "Strong"   ~ "Strong",
      partisanship == "Moderate" ~ "Moderate",
      TRUE ~ NA_character_
    )
  ) %>%
  filter(!is.na(condone))

prop_block <- function(df){
  take <- mean(df$condone[df$party_condition == "Inparty"],  na.rm = TRUE) * 100
  give <- mean(df$condone[df$party_condition == "Outparty"], na.rm = TRUE) * 100
  x <- c(sum(df$condone[df$party_condition=="Outparty"], na.rm = TRUE),
         sum(df$condone[df$party_condition=="Inparty"],  na.rm = TRUE))
  n <- c(sum(df$party_condition=="Outparty", na.rm = TRUE),
         sum(df$party_condition=="Inparty",  na.rm = TRUE))
  pval <- if (all(n > 0)) suppressWarnings(prop.test(x, n, alternative = "greater", correct = FALSE)$p.value) else NA_real_
  list(take = take, give = give, gap = give - take, pval = pval)
}

subs <- tibble::tribble(
  ~party_family, ~part_strength, ~row_lab,                 ~ord,
  "Republican",  "Strong",       "Strong Republicans",      1L,
  "Republican",  "Moderate",     "Moderate Republicans",    2L,
  "Republican",  NA_character_,  "Republicans (All)",       3L,
  "Democrat",    "Strong",       "Strong Democrats",        1L,
  "Democrat",    "Moderate",     "Moderate Democrats",      2L,
  "Democrat",    NA_character_,  "Democrats (All)",         3L
)

rows <- lapply(seq_len(nrow(subs)), function(i){
  sf  <- subs$party_family[i]
  st  <- subs$part_strength[i]
  d0  <- hp %>% filter(party_family == sf)
  d1  <- if (is.na(st)) d0 else d0 %>% filter(part_strength == st)
  s   <- prop_block(d1)
  tibble::tibble(
    party_family = sf,
    row_lab      = subs$row_lab[i],
    ord          = subs$ord[i],
    Take         = s$take,   # Inparty %
    Give         = s$give,   # Outparty %
    Gap          = s$gap,
    p_value      = s$pval
  )
})
fig2_data <- bind_rows(rows) %>%
  mutate(
    stars = case_when(
      is.na(p_value)  ~ "ns",
      p_value < 0.001 ~ "***",
      p_value < 0.01  ~ "**",
      p_value < 0.05  ~ "*",
      p_value < 0.10  ~ ".",
      TRUE            ~ "ns"
    ),
    Take_r = round(Take), Give_r = round(Give), Gap_r = round(Gap),
    GapLabel = paste0(Gap_r, " % ", stars),
    row_lab  = factor(row_lab,
                      levels = c("Strong Republicans","Moderate Republicans","Republicans (All)",
                                 "Strong Democrats","Moderate Democrats","Democrats (All)"))
  )

color_light_red   <- "#F2B8B5"
color_light_blue  <- "#B8C4E2"
color_red         <- "#E34A33"
color_blue        <- "#4F6DB8"
color_gray_axis   <- "#808080"

fig2_data <- fig2_data %>%
  mutate(
    barColor   = ifelse(party_family=="Democrat", color_light_blue, color_light_red),
    pointColor = ifelse(party_family=="Democrat", color_blue,       color_red)
  )

library(patchwork)
make_panel <- function(pfam, panel_title){
  dat <- fig2_data %>% filter(party_family == pfam) %>% arrange(ord)
  dat$y <- factor(dat$row_lab, levels = rev(levels(dat$row_lab)))
  
  p_main <- ggplot(dat) +
    geom_segment(aes(y = y, yend = y, x = Take, xend = Give, color = I(barColor)), size = 4) +
    geom_point(aes(y = y, x = Take, color = I(pointColor)), size = 5) +
    geom_point(aes(y = y, x = Give, color = I(pointColor)), size = 5) +
    geom_label(aes(y = y, x = Take, label = paste0("Inparty:", Take_r, "%")),
               vjust = -0.3, size = 3, label.size = NA, fill = NA, color = "black") +
    geom_label(aes(y = y, x = Give, label = paste0("Outparty:", Give_r, "%")),
               vjust = -0.3, size = 3, label.size = NA, fill = NA, color = "black") +
    scale_x_continuous(
      limits = c(5, 35),
      breaks = seq(5, 35, 5),
      labels = function(x) ifelse(x %% 10 == 0, paste0(x, "%"), "")
    ) +
    labs(title = panel_title, x = NULL, y = NULL) +
    theme_classic() +
    theme(
      plot.title        = element_text(face = "bold"),
      axis.text.x       = element_text(face = "bold", color = color_gray_axis),
      axis.text.y       = element_text(face = "bold"),
      legend.position   = "none",
      panel.border      = element_rect(fill = NA, color = "black", linewidth = 0.6)
    )
  
  p_gap <- ggplot(dat, aes(x = 0, y = y)) +
    geom_label(aes(label = GapLabel), fontface = "bold", size = 3.2,
               label.size = NA, fill = "white", alpha = 0.6) +
    theme_void() +
    ggtitle("Partisan gap") +
    theme(plot.title = element_text(face = "bold", size = 10, hjust = 0.5),
          panel.background = element_rect(fill = "white", colour = "white"),
          plot.margin = margin(t = 6, r = 6, b = 0, l = 0))
  
  p_main + p_gap + plot_layout(widths = c(5, 1))
}

p_rep <- make_panel("Republican", "Republicans")
p_dem <- make_panel("Democrat",   "Democrats")
Figure2 <- p_rep / p_dem + plot_layout(heights = c(1, 1))
Figure2

# Figure 3 - H3
Sum <- hypotheses_partisans %>%
  group_by(party_condition, disruption_severity) %>%
  summarise(
    mean = mean(response, na.rm=TRUE),
    se   = sd(response, na.rm=TRUE)/sqrt(n()),
    ci95 = se * qt(0.975, df = n()-1),
    ci83 = se * qt(0.917, df = n()-1),
    .groups = "drop"
  ) %>%
  mutate(
    party_condition = factor(party_condition, levels=c("No Party","Inparty","Outparty")),
    disruption_severity = factor(disruption_severity, levels=c("Mild","Severe"))
  )

Figure3 <- ggplot(Sum, aes(x = party_condition, y = mean, color = disruption_severity)) +
  geom_point(position = pd, size = 4) +
  geom_errorbar(aes(ymin = mean - ci95, ymax = mean + ci95),
                width = 0.0, position = pd, size = 0.8, alpha = 0.6) +
  geom_errorbar(aes(ymin = mean - ci83, ymax = mean + ci83),
                width = 0.0, position = pd, size = 2, alpha = 1) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        axis.title = element_text(face = "plain"),
        legend.position = "right",
        text = element_text(size = 11),
        axis.text.y = element_text(size = 11),
        axis.text.x = element_text(size = 12),
        plot.title = element_text(size = 12, face = "plain"),
        legend.text = element_text(size = 11)) +
  coord_cartesian(ylim = c(1.8, 3.6)) +
  labs(y = "Violence justification (full scale 1-7)", 
       x = "Victim Identity Marker", 
       color = "Provocation") +
  scale_color_manual(values = c("Mild" = "mediumorchid", "Severe" = "blueviolet"),
                     labels = c("Mild provocation", "Severe provocation"))
Figure3
ggsave("Figure_3_H3_means_severity.png", Figure3, width=7.5, height=4.5, dpi=300)
